Magnetostatics solution (linear case)ΒΆ
[1]:
%%capture
%run current_density.ipynb
[2]:
##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)
[3]:
R
[3]:
<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
with 23597 stored elements in Compressed Sparse Column format>
[7]:
order = 1
D = pde.int.assemble3(MESH, order = order)
phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)
K_Hcurl = curlphix_Hcurl @ D @ curlphix_Hcurl.T +\
curlphiy_Hcurl @ D @ curlphiy_Hcurl.T +\
curlphiz_Hcurl @ D @ curlphiz_Hcurl.T
KR = R.T@K_Hcurl@R
r = jx_L2 @ D @ phix_Hcurl.T +\
jy_L2 @ D @ phiy_Hcurl.T +\
jz_L2 @ D @ phiz_Hcurl.T
cholKR = chol(KR.tocsc())
A = R @ cholKR.solve_A(R.T@r)
curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
Bx = curlphix_Hcurl_P0.T @ A
By = curlphiy_Hcurl_P0.T @ A
Bz = curlphiz_Hcurl_P0.T @ A
[5]:
##############################################################################
# Storing to vtk
##############################################################################
grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')
[6]:
import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')
mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])
p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)
threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=1000)
p.add_mesh(arrows, color="orange")
mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")
p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")
[ ]: